1 Functions

library(rstatix)

2 Data

hmsc_signatue
 [1] "ANLN"         "TUBB6"        "AXL"          "IFT122"       "GFM1"         "SEPT2"        "PBRM1"        "INSIG1"       "CAPG"         "POP7"        
[11] "BHLHE40"      "MCM5"         "TRAPPC3"      "MRPS6"        "GINS2"        "NARS"         "LSM5"         "GUSB"         "EZR"          "FZD6"        
[21] "STAT2"        "TMPO"         "SLC35F5"      "MVB12A"       "FBLN1"        "COG4"         "MRPL20"       "GSDMD"        "VCL"          "MT1E"        
[31] "TUBA1C"       "ATXN10"       "USP39"        "MSH6"         "SDC2"         "OAT"          "SCAP"         "PRELID3B"     "KNTC1"        "WWTR1"       
[41] "ARF6"         "SMC5"         "DCTN2"        "AAMP"         "CSNK1A1"      "TBCB"         "KLF10"        "TPP1"         "EIF5A"        "IFT52"       
[51] "ETS2"         "EWSR1"        "NOP10"        "KRCC1"        "ASL"          "IGSF9"        "PHF5A"        "XRCC5"        "RNF168"       "UBE2T"       
[61] "IK"           "AC007250.3"   "YES1"         "SNRPD1"       "C9orf3"       "FDFT1"        "LASP1"        "SNRPC"        "PACSIN2"      "MTHFD2"      
[71] "PMVK"         "OBP2B"        "AKR1B1"       "SYPL1"        "LPAR2"        "PHB"          "THOC6"        "TMX1"         "SDHA"         "NUBP2"       
[81] "NDUFS7"       "TAOK3"        "MRPS9"        "PTRF"         "PSME2"        "MED15"        "PES1"         "RP11-19E11.1" "RP5-1102E8.2" "SEH1L"       
[91] "MCM3"         "POLR2L"       "DNAJC8"       "COPS6"        "PRKCDBP"      "KRT16"        "RNASEH2A"     "PDIA6"       
# get hpv genes and calc score
hpv_16_genes = rownames(opscc_hpvPos)[ startsWith(x =rownames(opscc_hpvPos),prefix =  "HPV16")]
hpv_score = FetchData(object = opscc_hpvPos,vars = hpv_16_genes,slot = "counts") %>% rowSums()
hpv_score_df = data.frame(hpv_score,row.names = names(hpv_score))
opscc_hpvPos = AddMetaData(object = opscc_hpvPos,metadata = hpv_score_df)

#set HPV to consider + if hpv_score>3
opscc_hpvPos$hpv = FetchData(object = opscc_hpvPos,vars = "hpv") %>%
 mutate(hpv = if_else(condition = hpv_score>=3,true = "HPV+",false = "HPV-")) %>% 
  pull(hpv)
opscc_hpvPos = ScaleData(object = opscc_hpvPos,features = hmsc_signatue)
Centering and scaling data matrix

  |                                                                                                                                                       
  |                                                                                                                                                 |   0%
  |                                                                                                                                                       
  |=================================================================================================================================================| 100%
DoHeatmap(object = opscc_hpvPos,features = hmsc_signatue,group.by = "hpv")

data = FetchData(object = opscc_hpvPos, vars = c(hmsc_signatue,"patient"))
data[] <- data %>%  group_by(patient) %>% mutate(across(is.numeric, ~ as.numeric(scale(.)))) %>% ungroup() 
data = data %>% select(-patient) 
data = data %>% t() %>% as.data.frame()

annotation_data = FetchData(object = opscc_hpvPos, vars = c("patient","hpv")) %>% dplyr::arrange(patient,hpv)
column_ha = HeatmapAnnotation(hpv = annotation_data[,2])
data = data[,rownames(annotation_data)] #order data like annotation
data = na.omit(data)
p = ComplexHeatmap::Heatmap(
  data,
  show_column_names = F,
  row_names_gp = grid::gpar(fontsize = 7),
  cluster_rows = T,
  top_annotation = column_ha,
  name = "Z-score expression",use_raster = F,cluster_columns = F,column_split   = annotation_data$patient,
  column_gap = unit(2, "mm"), border = TRUE,show_parent_dend_line = FALSE, show_column_dend = FALSE,cluster_column_slices = F
)
Warning: The input is a data frame, convert it to a matrix.The automatically generated colors map from the minus and plus 99^th of the absolute values in the matrix. There are outliers in the matrix
whose patterns might be hidden by this color mapping. You can manually set the color to `col` argument.

Use `suppressMessages()` to turn off this message.
  
p


# add score
hpv_signature_score = FetchData(object = opscc_hpvPos,vars = hmsc_signatue,slot = "data") %>% rowMeans()
Warning: The following requested variables were not found: AC007250.3, RP5-1102E8.2
hpv_signature = data.frame(hpv_signature_score,row.names = names(hpv_signature_score))
opscc_hpvPos = AddMetaData(object = opscc_hpvPos,metadata = hpv_signature)

genes_by_tp = FetchData(object = opscc_hpvPos,vars = c("hpv", "hpv_signature_score","patient"),slot = "data") 



#plot and split by patient:
stat.test <- genes_by_tp %>%
  group_by(patient) %>% 
  # filter (patient != "OP13") %>% 
  t_test(hpv_signature_score ~ hpv) %>%
  adjust_pvalue(method = "fdr") %>%
  add_significance() %>%
  add_xy_position(x = "hpv") %>%
  mutate(p.adj = signif(p.adj, digits = 2))
stat.test


plt = ggplot(genes_by_tp, aes(x = hpv, y = hpv_signature_score)) +
  geom_violin(scale = "width",aes(fill = hpv)) +
  geom_boxplot(width = 0.2, outlier.shape = NA) +
  stat_pvalue_manual(stat.test, label = "{p.adj}")  + #add p value
  facet_wrap("patient")+
  scale_y_continuous(expand = expansion(mult = c(0.05, 0.2))) + # Add 10% spaces between the p-value labels and the plot border
  ggtitle("HMSC HPV signature")
plt

LS0tCnRpdGxlOiAnYHIgcnN0dWRpb2FwaTo6Z2V0U291cmNlRWRpdG9yQ29udGV4dCgpJHBhdGggJT4lIGJhc2VuYW1lKCkgJT4lIGdzdWIocGF0dGVybiA9ICJcXC5SbWQiLHJlcGxhY2VtZW50ID0gIiIpYCcgCmF1dGhvcjogIkF2aXNoYWkgV2l6ZWwiCmRhdGU6ICdgciBTeXMudGltZSgpYCcKb3V0cHV0OiAKICBodG1sX25vdGVib29rOiAKICAgIGNvZGVfZm9sZGluZzogaGlkZQogICAgdG9jOiB5ZXMKICAgIHRvY19jb2xsYXBzZTogeWVzCiAgICB0b2NfZmxvYXQ6IAogICAgICBjb2xsYXBzZWQ6IEZBTFNFCiAgICBudW1iZXJfc2VjdGlvbnM6IHRydWUKICAgIHRvY19kZXB0aDogMQogICAgCiAgaHRtbF9kb2N1bWVudDogCiAgICBjb2RlX2ZvbGRpbmc6IGhpZGUKICAgIHRvYzogeWVzCiAgICB0b2NfY29sbGFwc2U6IHllcwogICAgdG9jX2Zsb2F0OiAKICAgICAgY29sbGFwc2VkOiBUUlVFCiAgICBudW1iZXJfc2VjdGlvbnM6IHRydWUKICAgIHRvY19kZXB0aDogMgogICAgZGYtcHJpbnQ6IHBhZ2VkCiAgICAKcGFyYW1zOgogIGRhdGFfb3V0X2RpcjogTlVMTAogIGZpZ3Nfb3V0X2RpcjogTlVMTAotLS0KCgoKIyBGdW5jdGlvbnMKCmBgYHtyIHdhcm5pbmc9RkFMU0V9CmxpYnJhcnkocnN0YXRpeCkKYGBgCgojIERhdGEKCmBgYHtyfQpvcHNjYyA9IHJlYWRSRFMoIi4vUmVuZGVyZWRfbm90ZWJvb2tzL0hQVl9PUFNDQ19BbmFseXNpc19QTUMxMDE5MTYzNC8wMV9wcmVwcm9jZXNzL29wc2NjX3N1ZXJhdC5SRFMiKQpsaWJyYXJ5KHJlYWR4bCkKbWV0YWRhdGEgPSByZWFkX2V4Y2VsKHBhdGggPSAiLi9JbnB1dF9kYXRhL0hQVl9PUFNDQ19BbmFseXNpc19QTUMxMDE5MTYzNC9OSUhNUzE4OTI3NTMtc3VwcGxlbWVudC1TdXBwX1RhYmxlcy54bHN4IiwKICAgICAgICAgIHNoZWV0ID0gIlM0IC0gQ2VsbCBUYWJsZSIscHJvZ3Jlc3MgPSBULHNraXAgPSAxLGNvbF9uYW1lcyA9IFQpCgojIHJlbW92ZSBIUFYtIHBhdGllbnRzCmhwdl9uZWc8LWMoIk9QMTAiLCJPUDEyIiwiT1AxNiIsIk9QMTkiLCJPUDgiKQpocHZfcG9zID0gdW5pcXVlKG1ldGFkYXRhJFBhdGllbnQpWyEgdW5pcXVlKG1ldGFkYXRhJFBhdGllbnQpICVpbiUgaHB2X25lZ10Kb3BzY2NfaHB2UG9zID0gc3Vic2V0KG9wc2NjLHN1YnNldCA9IHBhdGllbnQgJWluJSBocHZfcG9zKQoKaG1zY19kZWcgPSByZWFkLnRhYmxlKGZpbGUgPSAiLi9SZW5kZXJlZF9ub3RlYm9va3MvSE1TQy9IUFZfYW5hbHlzaXMvSFBWX2FuYWx5c2lzXzIwMjRfMTFfMTMvaHB2X2RlZ19kZi5jc3YiLHJvdy5uYW1lcyA9IDEpCmhtc2Nfc2lnbmF0dWUgPSByb3duYW1lcyhobXNjX2RlZ1tobXNjX2RlZyRmZHI8MC4wNSAmIGhtc2NfZGVnJGF2Z19sb2cyRkM+bG9nMigxLjMpLF0pCmhtc2Nfc2lnbmF0dWUKYGBgCgoKYGBge3J9CiMgZ2V0IGhwdiBnZW5lcyBhbmQgY2FsYyBzY29yZQpocHZfMTZfZ2VuZXMgPSByb3duYW1lcyhvcHNjY19ocHZQb3MpWyBzdGFydHNXaXRoKHggPXJvd25hbWVzKG9wc2NjX2hwdlBvcykscHJlZml4ID0gICJIUFYxNiIpXQpocHZfc2NvcmUgPSBGZXRjaERhdGEob2JqZWN0ID0gb3BzY2NfaHB2UG9zLHZhcnMgPSBocHZfMTZfZ2VuZXMsc2xvdCA9ICJjb3VudHMiKSAlPiUgcm93U3VtcygpCmhwdl9zY29yZV9kZiA9IGRhdGEuZnJhbWUoaHB2X3Njb3JlLHJvdy5uYW1lcyA9IG5hbWVzKGhwdl9zY29yZSkpCm9wc2NjX2hwdlBvcyA9IEFkZE1ldGFEYXRhKG9iamVjdCA9IG9wc2NjX2hwdlBvcyxtZXRhZGF0YSA9IGhwdl9zY29yZV9kZikKYGBgCgpgYGB7cn0KCiNzZXQgSFBWIHRvIGNvbnNpZGVyICsgaWYgaHB2X3Njb3JlPjMKb3BzY2NfaHB2UG9zJGhwdiA9IEZldGNoRGF0YShvYmplY3QgPSBvcHNjY19ocHZQb3MsdmFycyA9ICJocHYiKSAlPiUKIG11dGF0ZShocHYgPSBpZl9lbHNlKGNvbmRpdGlvbiA9IGhwdl9zY29yZT49Myx0cnVlID0gIkhQVisiLGZhbHNlID0gIkhQVi0iKSkgJT4lIAogIHB1bGwoaHB2KQoKYGBgCgoKYGBge3J9Cm9wc2NjX2hwdlBvcyA9IFNjYWxlRGF0YShvYmplY3QgPSBvcHNjY19ocHZQb3MsZmVhdHVyZXMgPSBobXNjX3NpZ25hdHVlKQpEb0hlYXRtYXAob2JqZWN0ID0gb3BzY2NfaHB2UG9zLGZlYXR1cmVzID0gaG1zY19zaWduYXR1ZSxncm91cC5ieSA9ICJocHYiKQpgYGAKYGBge3IgZmlnLmhlaWdodD0xNX0KZGF0YSA9IEZldGNoRGF0YShvYmplY3QgPSBvcHNjY19ocHZQb3MsIHZhcnMgPSBjKGhtc2Nfc2lnbmF0dWUsInBhdGllbnQiKSkKZGF0YVtdIDwtIGRhdGEgJT4lICBncm91cF9ieShwYXRpZW50KSAlPiUgbXV0YXRlKGFjcm9zcyhpcy5udW1lcmljLCB+IGFzLm51bWVyaWMoc2NhbGUoLikpKSkgJT4lIHVuZ3JvdXAoKSAKZGF0YSA9IGRhdGEgJT4lIHNlbGVjdCgtcGF0aWVudCkgCmRhdGEgPSBkYXRhICU+JSB0KCkgJT4lIGFzLmRhdGEuZnJhbWUoKQoKYW5ub3RhdGlvbl9kYXRhID0gRmV0Y2hEYXRhKG9iamVjdCA9IG9wc2NjX2hwdlBvcywgdmFycyA9IGMoInBhdGllbnQiLCJocHYiKSkgJT4lIGRwbHlyOjphcnJhbmdlKHBhdGllbnQsaHB2KQpjb2x1bW5faGEgPSBIZWF0bWFwQW5ub3RhdGlvbihocHYgPSBhbm5vdGF0aW9uX2RhdGFbLDJdKQpkYXRhID0gZGF0YVsscm93bmFtZXMoYW5ub3RhdGlvbl9kYXRhKV0gI29yZGVyIGRhdGEgbGlrZSBhbm5vdGF0aW9uCmRhdGEgPSBuYS5vbWl0KGRhdGEpCnAgPSBDb21wbGV4SGVhdG1hcDo6SGVhdG1hcCgKICBkYXRhLAogIHNob3dfY29sdW1uX25hbWVzID0gRiwKICByb3dfbmFtZXNfZ3AgPSBncmlkOjpncGFyKGZvbnRzaXplID0gNyksCiAgY2x1c3Rlcl9yb3dzID0gVCwKICB0b3BfYW5ub3RhdGlvbiA9IGNvbHVtbl9oYSwKICBuYW1lID0gIlotc2NvcmUgZXhwcmVzc2lvbiIsdXNlX3Jhc3RlciA9IEYsY2x1c3Rlcl9jb2x1bW5zID0gRixjb2x1bW5fc3BsaXQgICA9IGFubm90YXRpb25fZGF0YSRwYXRpZW50LAogIGNvbHVtbl9nYXAgPSB1bml0KDIsICJtbSIpLCBib3JkZXIgPSBUUlVFLHNob3dfcGFyZW50X2RlbmRfbGluZSA9IEZBTFNFLCBzaG93X2NvbHVtbl9kZW5kID0gRkFMU0UsY2x1c3Rlcl9jb2x1bW5fc2xpY2VzID0gRgopCiAgCnAKYGBgCmBgYHtyfQoKIyBhZGQgc2NvcmUKaHB2X3NpZ25hdHVyZV9zY29yZSA9IEZldGNoRGF0YShvYmplY3QgPSBvcHNjY19ocHZQb3MsdmFycyA9IGhtc2Nfc2lnbmF0dWUsc2xvdCA9ICJkYXRhIikgJT4lIHJvd01lYW5zKCkKaHB2X3NpZ25hdHVyZSA9IGRhdGEuZnJhbWUoaHB2X3NpZ25hdHVyZV9zY29yZSxyb3cubmFtZXMgPSBuYW1lcyhocHZfc2lnbmF0dXJlX3Njb3JlKSkKb3BzY2NfaHB2UG9zID0gQWRkTWV0YURhdGEob2JqZWN0ID0gb3BzY2NfaHB2UG9zLG1ldGFkYXRhID0gaHB2X3NpZ25hdHVyZSkKCmBgYAoKCmBgYHtyIGZpZy5oZWlnaHQ9NiwgZmlnLndpZHRoPTh9CgpnZW5lc19ieV90cCA9IEZldGNoRGF0YShvYmplY3QgPSBvcHNjY19ocHZQb3MsdmFycyA9IGMoImhwdiIsICJocHZfc2lnbmF0dXJlX3Njb3JlIiwicGF0aWVudCIpLHNsb3QgPSAiZGF0YSIpIAoKCgojcGxvdCBhbmQgc3BsaXQgYnkgcGF0aWVudDoKc3RhdC50ZXN0IDwtIGdlbmVzX2J5X3RwICU+JQogIGdyb3VwX2J5KHBhdGllbnQpICU+JSAKICAjIGZpbHRlciAocGF0aWVudCAhPSAiT1AxMyIpICU+JSAKICB0X3Rlc3QoaHB2X3NpZ25hdHVyZV9zY29yZSB+IGhwdikgJT4lCiAgYWRqdXN0X3B2YWx1ZShtZXRob2QgPSAiZmRyIikgJT4lCiAgYWRkX3NpZ25pZmljYW5jZSgpICU+JQogIGFkZF94eV9wb3NpdGlvbih4ID0gImhwdiIpICU+JQogIG11dGF0ZShwLmFkaiA9IHNpZ25pZihwLmFkaiwgZGlnaXRzID0gMikpCnN0YXQudGVzdAoKCnBsdCA9IGdncGxvdChnZW5lc19ieV90cCwgYWVzKHggPSBocHYsIHkgPSBocHZfc2lnbmF0dXJlX3Njb3JlKSkgKwogIGdlb21fdmlvbGluKHNjYWxlID0gIndpZHRoIixhZXMoZmlsbCA9IGhwdikpICsKICBnZW9tX2JveHBsb3Qod2lkdGggPSAwLjIsIG91dGxpZXIuc2hhcGUgPSBOQSkgKwogIHN0YXRfcHZhbHVlX21hbnVhbChzdGF0LnRlc3QsIGxhYmVsID0gIntwLmFkan0iKSAgKyAjYWRkIHAgdmFsdWUKICBmYWNldF93cmFwKCJwYXRpZW50IikrCiAgc2NhbGVfeV9jb250aW51b3VzKGV4cGFuZCA9IGV4cGFuc2lvbihtdWx0ID0gYygwLjA1LCAwLjIpKSkgKyAjIEFkZCAxMCUgc3BhY2VzIGJldHdlZW4gdGhlIHAtdmFsdWUgbGFiZWxzIGFuZCB0aGUgcGxvdCBib3JkZXIKICBnZ3RpdGxlKCJITVNDIEhQViBzaWduYXR1cmUiKQpwbHQKYGBgCg==